Description

This script merges CalMAPPER activity data with treatment polygons. This is a necessary step before analyzing prescribed fire data in CalMAPPER.

Define path to Calmapper and to the gdb file

calmapper_dir <- fs::dir_ls(ref_path, recurse = T, glob = '*CalMAPPER', type = 'directory')
calmapper_gdb <- fs::dir_ls(calmapper_dir, recurse = T, glob = '*.gdb')
if(length(calmapper_gdb) > 1){
  stop( "There's more than one CalMAPPER .gdb file. script assumes only 1.")
}
st_layers(calmapper_gdb)
## Driver: OpenFileGDB 
## Available layers:
##                 layer_name     geometry_type features fields                 crs_name
## 1 CMDash_ProjectTreatments     Multi Polygon     1560     21 WGS 84 / Pseudo-Mercator
## 2     CMDash_TreatmentPols     Multi Polygon     2784     23 WGS 84 / Pseudo-Mercator
## 3    CMDash_TreatmentLines Multi Line String       23     22 WGS 84 / Pseudo-Mercator
## 4     CMDash_TreatmentPnts       Multi Point      136     20 WGS 84 / Pseudo-Mercator
## 5        CMDash_Activities                NA    20439     38                     <NA>
## 6          CMDash_Metadata                NA        1      7                     <NA>
## 7       CalMapper_fire_act     Multi Polygon     2799     39 WGS 84 / Pseudo-Mercator
act <- st_read(calmapper_gdb, layer = 'CMDash_Activities')
## Reading layer `CMDash_Activities' from data source 
##   `C:\Users\ctubbesi\OneDrive - California Air Resources Board\Documents\Reference data\CALFIRE spatial data\CalMAPPER\CALFIRE_FuelReductionProjects_2023\CALFIRE_FuelReductionProjects.gdb' 
##   using driver `OpenFileGDB'
## Warning: no simple feature geometries present: returning a data.frame or tbl_df
trt <- st_read(calmapper_gdb, layer = 'CMDash_TreatmentPols')
## Reading layer `CMDash_TreatmentPols' from data source 
##   `C:\Users\ctubbesi\OneDrive - California Air Resources Board\Documents\Reference data\CALFIRE spatial data\CalMAPPER\CALFIRE_FuelReductionProjects_2023\CALFIRE_FuelReductionProjects.gdb' 
##   using driver `OpenFileGDB'
## Simple feature collection with 2784 features and 23 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13842530 ymin: 3842624 xmax: -12941530 ymax: 5161286
## Projected CRS: WGS 84 / Pseudo-Mercator

Clean

Remove non-treatment “activities”

act <- act %>% 
  filter(!ACTIVITY_DESCRIPTION %in% c("GIS Validation", "Education Outreach", "Project Administration", "Planning Meeting", "Public Contacts", "Water Site Development"))

Merge Rx fire Activity data with treatment polygons

act_sf <- right_join(trt %>% select(PROJECT_ID, TREATMENT_ID), 
           act)
## Joining with `by = join_by(PROJECT_ID, TREATMENT_ID)`

Filter to fire

Find relevant treatment names

act_sf %>% 
  st_drop_geometry() %>% 
  group_by(ACTIVITY_DESCRIPTION) %>% 
  count() %>% 
  print(n=50)
## # A tibble: 47 × 2
## # Groups:   ACTIVITY_DESCRIPTION [47]
##    ACTIVITY_DESCRIPTION                      n
##    <chr>                                 <int>
##  1 Air Curtain Burner                       13
##  2 Biomass Removal (Bone Dry Tons)          85
##  3 Boundary Mapping                         43
##  4 Broadcast Burn                          868
##  5 Chaining                                100
##  6 Chipping                               2216
##  7 Commercial Thinning (Cable Yarding)      20
##  8 Commercial Thinning (Tractor Yarding)    42
##  9 Crushing                                 51
## 10 Cultural Burning                          5
## 11 Dozer Line                               35
## 12 Environmental Review                    110
## 13 Erosion Control                          12
## 14 Follow up - Herbicide                    58
## 15 Follow up - Other                         6
## 16 Follow up - Slash disposal              177
## 17 Fuel Break (Shaded)                      55
## 18 Grazing                                  52
## 19 Hand Line                               108
## 20 Herbicide (Post-Treatment)               43
## 21 Herbicide (Pre-Treatment)                17
## 22 Invasive Plant Removal                    7
## 23 Land Conservation                         4
## 24 Limbing and Bucking                    1204
## 25 Lop and Scatter                         717
## 26 Mastication                            1243
## 27 Pile Burning                           1926
## 28 Piling (Manual)                        2533
## 29 Piling (Mechanical)                     375
## 30 Pruning                                 376
## 31 Public Meetings                          21
## 32 RPF Supervision                          86
## 33 Rangeland Mowing                         61
## 34 Release - Herbicide                      16
## 35 Release - Mechanical                     99
## 36 Release - Other                          10
## 37 Road Grading                              6
## 38 Site Assessment                          93
## 39 Site Preparation (CFIP)                  71
## 40 Site Preparation (Manual)                77
## 41 Site Preparation (Mechanical)             7
## 42 Site Preparation (RxBurn)               187
## 43 Thinning                                287
## 44 Thinning (Manual)                      4400
## 45 Thinning (Mechanical)                   442
## 46 Trees Felled (> 6in dbh)               1026
## 47 <NA>                                      1

Filter by treatment

act_sf_fire <- 
  act_sf %>% 
  filter(ACTIVITY_DESCRIPTION %in% c("Broadcast Burn", "Cultural Burning", "Pile Burning"))

Create DURATION and YEAR columns

DURATION

act_sf_fire <- act_sf_fire %>% 
  mutate(DURATION = as.Date(ACTIVITY_END) - as.Date(ACTIVITY_START)+1)

YEAR

act_sf_fire <- act_sf_fire %>% 
  mutate(YEAR = year(ACTIVITY_END))

Check

act_sf_fire %>% 
  select(ACTIVITY_START, ACTIVITY_END, DURATION, YEAR) %>% 
  head() %>% 
  nrow()
## [1] 6

Save

save(act_sf_fire, file = "Rdata/CalMapper_activities_fire.Rdata")
write.csv(act_sf_fire %>% st_drop_geometry(), file = "~/Reference data/CALFIRE spatial data/CalMapper/activities_fire.csv", row.names = F)
write.csv(trt %>% st_drop_geometry(), file = "~/Reference data/CALFIRE spatial data/CalMapper/treatments_fire.csv", row.names = F)

Are there Rx fires recorded in Treatments data but not in Activities data or vice versa?

act_broadcast <- act_sf_fire %>% 
  filter(ACTIVITY_DESCRIPTION == "Broadcast Burn")
trt %>% 
  st_drop_geometry() %>% 
  group_by(TREATMENT_OBJECTIVE) %>% 
  count()
## # A tibble: 5 × 2
## # Groups:   TREATMENT_OBJECTIVE [5]
##   TREATMENT_OBJECTIVE        n
##   <chr>                  <int>
## 1 Broadcast Burn           559
## 2 Forestland Stewardship   296
## 3 Fuel Break               215
## 4 Fuel Reduction          1542
## 5 Right of Way Clearance   172
trt_broadcast_complete <- trt %>% 
  filter(TREATMENT_OBJECTIVE == "Broadcast Burn") %>% 
  filter(ACTIVITY_STATUS == "Complete")

Check if there are treatment IDs in act_broadcast that aren’t in trt – this shows whether the merge was complete

act_broadcast %>% 
  filter(!TREATMENT_ID %in% trt$TREATMENT_ID) %>% 
  nrow()
## [1] 0
act_broadcast %>% 
  filter(!TREATMENT_NAME %in% trt$TREATMENT_NAME) %>% 
  nrow()
## [1] 0

Check if there are treatment IDs in completed broadcast treatments that aren’t in activities

trt_broadcast_complete %>% 
  filter(!TREATMENT_ID %in% act_broadcast$TREATMENT_ID) %>% 
  mapview()
trt_broadcast_complete %>% 
  filter(!TREATMENT_NAME %in% act_broadcast$TREATMENT_NAME) %>% 
  nrow()
## [1] 7

Export those to put in .ppt

missing_from_activities <- trt_broadcast_complete %>% 
  filter(!TREATMENT_ID %in% act_broadcast$TREATMENT_ID) %>% 
  select(CALMAPPER_ID, PROJECT_NAME, TREATMENT_NAME, TREATMENT_OBJECTIVE, PROJECT_STATUS, ACTIVITY_STATUS, PROJECT_START_DATE, PROJECT_END_DATE, TREATMENTAREA_ACRES) %>% 
  st_drop_geometry()
missing_from_activities %>% 
  summarize(missing_acres = sum(TREATMENTAREA_ACRES))
##   missing_acres
## 1        653.99
write_excel_csv(missing_from_activities, "calmapper_trt_not_act.xls")

Visualize broadcast treatments with no overlapping activities

map <- mapview(list(trt_broadcast_complete, act_broadcast), col.regions=list("red","blue"),col=list("red","blue"))
## Warning in clean_columns(as.data.frame(obj), factorsAsCharacter): Dropping column(s) DURATION of class(es) difftime
map

It’s safe to use the Activities information and not the Treatment information.

Write activities data to gdb

st_write(act_sf_fire, dsn = calmapper_gdb, layer = "CalMapper_fire_act", delete_layer = T)
## Warning in clean_columns(as.data.frame(obj), factorsAsCharacter): Dropping column(s) DURATION of class(es) difftime
## Deleting layer `CalMapper_fire_act' using driver `OpenFileGDB'
## Writing layer `CalMapper_fire_act' to data source 
##   `C:/Users/ctubbesi/OneDrive - California Air Resources Board/Documents/Reference data/CALFIRE spatial data/CalMAPPER/CALFIRE_FuelReductionProjects_2023/CALFIRE_FuelReductionProjects.gdb' using driver `OpenFileGDB'
## Writing 2799 features with 39 fields and geometry type Multi Polygon.

Remove duplicates from activities data

act_sf_fire %>% 
  filter(duplicated(SHAPE) & duplicated(YEAR) & duplicated(ACTIVITY_DESCRIPTION)) %>% 
  arrange(SHAPE) %>% 
  select(PROJECT_ID, TREATMENT_ID, ACTIVITY_ID, AQ_ID, PROJECT_NAME, TREATMENT_NAME, ACTIVITY_NAME, ACTIVITY_DESCRIPTION, )
## Simple feature collection with 1666 features and 8 fields (with 6 geometries empty)
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13836290 ymin: 3843456 xmax: -12941590 ymax: 5161029
## Projected CRS: WGS 84 / Pseudo-Mercator
## First 10 features:
##    PROJECT_ID TREATMENT_ID ACTIVITY_ID  AQ_ID                    PROJECT_NAME               TREATMENT_NAME     ACTIVITY_NAME
## 1         399          722       17318  73987 Mount Konocti Interface Project Clearlake Riviera Fuel Break 2020 Pile Burning
## 2         399          722       17318  80362 Mount Konocti Interface Project Clearlake Riviera Fuel Break 2020 Pile Burning
## 3         406          726       17104  64347   Las Posadas Shaded Fuel Break         Las Posadas Thinning 2019 Pile Burning
## 4         406          726       17314  67257   Las Posadas Shaded Fuel Break         Las Posadas Thinning 2020 Pile Burning
## 5         406          726       17314  81297   Las Posadas Shaded Fuel Break         Las Posadas Thinning 2020 Pile Burning
## 6         406          726       17314  74024   Las Posadas Shaded Fuel Break         Las Posadas Thinning 2020 Pile Burning
## 7        2023         3038       17326  78121                      Lake Curry               LAKE CURRY DAM 2020 Pile Burning
## 8        2023         3038       17326  95826                      Lake Curry               LAKE CURRY DAM 2020 Pile Burning
## 9        2023         3038       21708 173142                      Lake Curry               LAKE CURRY DAM 2021 Pile Burning
## 10       2023         3038       21708 173143                      Lake Curry               LAKE CURRY DAM 2021 Pile Burning
##    ACTIVITY_DESCRIPTION                          SHAPE
## 1          Pile Burning MULTIPOLYGON (((-13662768 4...
## 2          Pile Burning MULTIPOLYGON (((-13662768 4...
## 3          Pile Burning MULTIPOLYGON (((-13626585 4...
## 4          Pile Burning MULTIPOLYGON (((-13626585 4...
## 5          Pile Burning MULTIPOLYGON (((-13626585 4...
## 6          Pile Burning MULTIPOLYGON (((-13626585 4...
## 7          Pile Burning MULTIPOLYGON (((-13595922 4...
## 8          Pile Burning MULTIPOLYGON (((-13595922 4...
## 9          Pile Burning MULTIPOLYGON (((-13595922 4...
## 10         Pile Burning MULTIPOLYGON (((-13595922 4...

Save 2020 to use for Venn Diagram

act_2020 <- act_sf_fire %>% 
  filter(YEAR == 2020)

Mask out federal lands

dsn = paste(ref_path, "/National Forest Boundaries", sep = "")
NF <- st_read(dsn = dsn, "NF_CA_diss")
## Reading layer `NF_CA_diss' from data source 
##   `C:\Users\ctubbesi\OneDrive - California Air Resources Board\Documents\Reference data\National Forest Boundaries' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.0956 ymin: 32.65371 xmax: -116.3189 ymax: 42.00876
## Geodetic CRS:  NAD83
NF <- st_transform(NF, st_crs(act_2020))
save(NF,  file = "Rdata/NF.Rdata")
act_2020 <- st_difference(act_2020, NF)
## Warning: attribute variables are assumed to be spatially constant throughout all geometries

Save

#st_write(act_2020, dsn = "shapefiles", layer = "CalMapper_activities_fire_2020", driver = "ESRI Shapefile", delete_layer = T)
save(act_2020, file = "Rdata/CalMapper_activities_fire_2020.Rdata")